#!/usr/bin/perl -w
use strict;
use DBI;
use DBD::Pg;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use chromosomes;
use linewrap;

use db_parameters;
use ortho_species;

print <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    analyse_count_homo_regions_without_orthologs


	depend on
		orthologs.transcript_clusters
		homo.gene_loci
		canis.gene_loci


888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE

#
#	Connect to genome database
#
my $dbh = connect_to_panda();
my @db_param = get_db_parameters();

sub get_gaps($$\%\@\%)
{
	my ($dbh, $species, $gap_cnts, $gap_sizes, $gap_gene_ids) = @_;
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
        DISTINCT (GENE_ID), 
        a.gene_id is null as missing, 
        chromosome, 
        gene_position, 
        start,
        finish 
    FROM 
        (
            SELECT
                gene_id
            FROM 
                orthologs.ortholog_sets
            WHERE
                orthologs.ortholog_sets.protocol_id = 1 AND
                NOT ortholog_type ~ 'orphan'
        ) AS a RIGHT JOIN 
        $species.gene_loci USING (gene_id)
    ORDER BY gene_position;
PL/SQLCMD

	my $prev_chrm = '';
	my $in_gap = 0;
	my $gap_cnt = 0;
	my $gap_beg = 0;
	my $gap_end = 0;
	my $gap_beg_line = 0;
	my $gap_end_line = 0;
	my $gap_beg_gene = '';
	my $gap_end_gene = '';
	my $gap_beg_gp = 0;
	my $gap_end_gp = 0;

	my $data = $dbh->selectall_arrayref($sql_cmd);
	#for (@$data)
	#{
	#	print STDERR join ("\t",  @$_),  "\n";
	#}
	my @gap_gene_ids;
	for my $gene(@$data)
	{
		my ($gene_id, $missing, $chrm, $gp, $start, $finish) = @$gene;
		if ($prev_chrm ne $chrm)
		{
			# close previous gap
			if ($in_gap)
			{
				push(@{$gap_gene_ids->{$gap_cnt}}, [@gap_gene_ids]);
				@gap_gene_ids = ();
				$gap_cnts->{$gap_cnt}++;
				die "$gap_beg > $gap_end close gap at end of chromosome" if $gap_beg > $gap_end;
				push(@$gap_sizes, $gap_end - $gap_beg);
			}

			$prev_chrm = $chrm;
			$in_gap = 0;
			$gap_cnt = 0;
			$gap_end = 0;
			$gap_beg = 0;
			$gap_beg_gene = $gene_id;
			$gap_beg_gp = $gp;
			$gap_beg_line = join "\t", @$gene;
			#print STDERR "New chromosome\n";
		}

		# extend gap
		if ($missing)
		{
			$in_gap = 1;
			#print STDERR "Missing\n";
			$gap_cnt++;
			# save gene in gap
			push(@gap_gene_ids, $gene_id);
			# gap ends either with last gene on chromosome (kludge) or with first gene after gap
			$gap_end = $finish;
			$gap_end_gene = $gene_id;
			$gap_end_gp = $gp;
			$gap_end_line = join "\t", @$gene;

		}

		#
		# close gap
		#
		if ($in_gap && !$missing)
		{
			#print STDERR "Close gap\n";
			# gap ends with first gene after gap
			$gap_end = $start;
			$gap_end_gene = $gene_id;
			$gap_end_gp = $gp;
			$gap_end_line = join "\t", @$gene;

			push(@{$gap_gene_ids->{$gap_cnt}}, [@gap_gene_ids]);
			@gap_gene_ids = ();
			$gap_cnts->{$gap_cnt}++;
			#die "$gap_beg > $gap_end close gap, $gap_beg_gp $gap_beg_gene to $gap_end_gp $gap_end_gene\n$gap_beg_line\n$gap_end_line\n" if $gap_beg > $gap_end;
			push(@$gap_sizes, $gap_end - $gap_beg) if $gap_end > $gap_beg;

			$gap_cnt = 0;

		}
		
		if (!$missing)
		{
			#print STDERR "Not missing\n";
			# if next gene is a gap then this finish starts the next gap
			$gap_beg = $finish;
			$gap_beg_gene = $gene_id;
			$gap_beg_gp = $gp;
			$gap_beg_line = join "\t", @$gene;
			$in_gap = 0;
		}



	}
}

sub print_genes_in_gaps(\%$$)
{
	my ($gap_gene_ids, $file_name, $min_gap_size) = @_;

	# make sure there are genes in large gaps
	my @gap_sizes = sort {$a <=> $b} keys %$gap_gene_ids;
	if ($gap_sizes[-1] < $min_gap_size)
	{
		print "\t\tNo genes in gaps larger than $min_gap_size\n";
		return;
	}

	# print lists of gene_ids in gaps
	my $cnt_gaps = 0;
	my $cnt_gap_genes = 0;
	open GAPGENES, ">$FindBin::Bin/$file_name" or die;

	for my $size(@gap_sizes)
	{
		next unless ($size >= $min_gap_size);
		for my $genes(@{$gap_gene_ids->{$size}})
		{
			$cnt_gaps++;
			$cnt_gap_genes +=@$genes;
			print GAPGENES  join ("\n",  @$genes), "\n";
		}
	}
	my $gap_str = $min_gap_size ?  ">= $min_gap_size genes in size" : "";
	print "\t\t$cnt_gaps\tgaps $gap_str\n";
	print "\t\t$cnt_gap_genes\tgenes in gaps $gap_str\n";
}


for my $species($ortho_name1, $ortho_name2)
{
	print "\tAnalyse gaps in $species...\n";

	open GAPS, ">$FindBin::Bin/$species.gap" or die;
	my (%gap_cnts, @gap_sizes, %gap_gene_ids);
	get_gaps($dbh, $species, %gap_cnts, @gap_sizes, %gap_gene_ids);
	print GAPS "size\tcount\n";
	for my $size(sort {$a <=> $b} keys %gap_cnts)
	{
		print GAPS "$size\t$gap_cnts{$size}\n";
	}
	print GAPS "\n";
	print GAPS join ("\n", sort {$a<=> $b} @gap_sizes);

	print "\t$species genes from gaps of >= 10\n";
	print_genes_in_gaps(%gap_gene_ids, "$species.biggap.genes", 10);
	print "\tAll $species genes from gaps\n";
	print_genes_in_gaps(%gap_gene_ids, "$species.gap.genes", 0);
	print "\n\n";
}

